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Abstract 

Order-preserving couplings are elegant tools for obtaining robust es- 
timates of the time-dependent and stationary distributions of Markov 
processes that are too complex to be analyzed exactly. The start- 
ing point of this paper is to study stochastic relations, which may be 
viewed as natural generalizations of stochastic orders. This generaliza- 
tion is motivated by the observation that for the stochastic ordering of 
two Markov processes, it suffices that the generators of the processes 
preserve some, not necessarily reflexive or transitive, sub-relation of the 
order relation. The main contributions of the paper are an algorith- 
mic characterization of stochastic relations between finite spaces, and 
a truncation approach for comparing infinite-state Markov processes. 
The methods are illustrated with applications to loss networks and 
parallel queues. 
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1 Introduction 



Comparison techniques based on stochastic orders [12, 13, 14] are key to 
obtaining upper and lower bounds for complicated random variables and 
processes in terms of simpler random elements. Consider for example two 
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ergodic discrete-time Markov processes X and Y with stationary distribu- 
tions fix and fj,y-> taking values in a common ordered state space, and denote 
by <st the corresponding stochastic order. Then the upper bound 

fJ-x < s t (1) 

can be established [G] without explicit knowledge of by verifying that 
the corresponding transition probability kernels Px and Py satisfy 

x<y Px(x,-) <st Py(y,-). (2) 

Analogous conditions for continuous-time Markov processes on countable 
spaces have been derived by Whitt [15] and Massey [11], and later extended 
to more general jump processes by Brandt and Last [2]. 

Less stringent sufficient conditions for obtaining (1) have recently been 
found using a new theory stochastic relations [9]. Two random variables 
are stochastically related, denoted by X ~ st Y, if there exists a coupling 
(X ,Y) of X and Y such that X ~ Y almost surely, where ~ denotes some 
relation between the state spaces of X and Y. The main motivation for 
this definition is that (2) is by no means necessary for (1); a less stringent 
sufficient condition is that 

x~y PxM ~BtiV(j/,-) (3) 

for some, not necessarily symmetric or transitive, nontrivial subrelation of 
the underlying order relation. Another advantage of the generalized defi- 
nition is that X and Y are no longer required to take values in the same 
state space, leading to greater flexibility in the search for bounding random 
elements Y . For example, to study whether f(X) < st g(Y) for some given 
real functions / and g defined on the state spaces of X and Y, we may define 
a relation x ~ y by the condition f(x) < g(y) [3]. 

The rest of the paper is outlined as follows. After recalling the basic 
definitions, Section 2 presents a numerical algorithm for verifying stochastic 
relations between finite spaces, together with an analysis of computational 
complexity. Section 3 recalls how a recursive subrelation algorithm may be 
used to find Markovian couplings preserving a given relation. In Section 4, 
a new truncation approach is presented that allows to precisely compute 
truncated outcomes of the subrelation algorithm for infinite-state Markov 
processes with locally bounded jumps. Section 5 discusses applications to 
loss networks and parallel queues, and Section 6 concludes the paper. 
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2 Stochastic relations 



2.1 General definitions 

We shall here recall the definitions of stochastic relations between countable 
spaces. Probability measures \i on a countable state space S shall be viewed 
as a probability vectors via identifying fj,(x) = fj,({x}). For a treatment on 
more general spaces, see [9]. 

A relation between Si and S2 is subset of Si x S2. Given a nontrivial 
(R 7^ 0) relation R between Si and S2, we write 

x ~ y, 

if (x, y) G i?. The relation may equivalently be viewed as a matrix so that 
i?(x, y) = 1 if x ~ y and i?(x, y) = otherwise. A coupling of probability 
vectors \i on Si and v on S2 is a probability vector A on S\ x S2 with 
marginals /i and v, that is, 



A(x,y) = /x(x) for all 1 G Si, 



y&S 2 

A(x,y) = i/(y) for all y 6 S2. 

For probability vectors fj, on Si and v on S2 we denote 

M ~st f, 

and say that is stochastically related to u, if there exists a coupling A of 
fi and f such that 

E M^y) = L 

The relation R s t = {(//, i>) : /x ~ s t i>} is called the stochastic relation gener- 
ated by R. Observe that two Dirac masses satisfy 5 X ~ s t o~y if and only if 
x ~ y. In this way the stochastic relation R st may be regarded as a natural 
randomization of the underlying relation R. 

The following result in [9], which is rephrased here for ease of refer- 
ence, provides an analytical method to check whether a pair of probability 
measures are stochastically related. 

Theorem 1. [9] Two probability vectors \x and v are stochastically related 
with respect to R if and only if 



max R(x,y) j 



<v) (4) 



xGU y£S2 
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for all finite U C S\ , or equivalently, if and only if 

waxR(x,y))ii(x) (5) 

y&V xGSi ^ V ' 

for all finite V C S2 ■ 

A random variable X is stochastically related to a random variable V, 
denoted by X ~ st V , if the distribution of X is stochastically related to the 
distribution of Y . Observe that X and Y do not need to be defined on the 
same probability space. Recall that a coupling of random variables X and Y 
is a bivariate random variable whose distribution couples the distributions 
of X and Y . Hence X ~ st Y if and only if there exists a coupling (X, Y) of 
X and Y such that X ~ Y" almost surely. 

Example 1. If < is an order (reflexive and transitive) relation on a space 
S, then the corresponding stochastic relation < st is called a stochastic order. 
Using Strassen's classical theorem [12], we see that X < st Y if and only if 
E/(X) < E/(Y") for all positive increasing functions / on S. 



2.2 Stochastic relations between finite spaces 

Let R be a relation between finite spaces S\ and £2, and denote by R s t the 
corresponding stochastic relation. Then Theorem 1 may be used to deter- 
mine whether [i ~ st v. However, this requires to check the inequality (4) 
for all subsets of Si, which is computationally infeasible unless the spaces 
are small. The following result shows that less checks may be sufficient. We 
shall denote the support of a probability vector fj, by Up = {x : fi(x) > 0}. 
Moreover, we denote by -F([/, Z2) the set of vectors with components in 
{0, 1} indexed by elements of U, which may also be identified as the set of 
all subsets of U. 

Theorem 2. Two probability vectors fi and v with supports and U v are 
stochastically related with respect to R if and only if 

V f(x)n(x) < V msx.[f(x)R(x,y)]u{y) (6) 

* — * — xGU^ 

for all f £ FfUfj,,^), or equivalently, if and only if 

y^5 , (y)^(y)< Y] max [R{x, y)g(y)] n(x) (7) 
3/ef/i/ xec/ M 

/or a// 5 G F(U V ,Z 2 ). 
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Proof. In light of Theorem 1, it suffices to show the equivalence of (4) 
and (6), and the equivalence of (5) and (7). Observe that (4) directly im- 
plies (6), because the members of F^U^'L-i) may be identified with the 
indicator functions of subsets of U^. To prove the converse, assume that (6) 
holds, and let U be an arbitrary subset of Si. Define f(x) = l(x G UCl U^). 
Then 



Hence (4) holds. Proving the equivalence of (5) and (7) is completely anal- 



Algorithm 1 Determining whether fi ~ s t v. 
U„ <- {x i Si : > 0} 
U v ^{y€S 2 : u(y) > 0} 
if > #U V then 

flip fx <-> v, Ufj, <-> C/jy, Si <-> S 2 
end if 
6 <— true 

for i = l,...,2* u » do 

/ <— i-th vector in FiTJ^^TL?) 

u r <- Ej/Gf/, [ max xec/ M /(x)i?(x,y)] 
if ^ > f r then 
b <— false 
break 
end if 
end for 
return b 



Algorithm 1 describes how Theorem 2 can be applied to numerically de- 
termine whether \i ~ st v. The interchange of the variables in the beginning 




ogous. 



□ 
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corresponds to using (6) if the support of (i is smaller than is, and (7) oth- 
erwise. Inspection of Algorithm 1 shows that the computational complexity 
of determining whether fi ~ s t v is of the order 



where n'i and n' 2 denote the cardinalities of the supports of /i and v. The 
algorithm is very slow when both state spaces are large and fi and v have 
positive mass in all states. However, in many applications related to struc- 
tured Markov chains we may assume that the \x and v have small supports 
(see Section 5). 

Remark 1. The verification of [i ~ s t v can be carried out faster, if the un- 
derlying relation R has some structure that can be employed. For example, 
for the natural order on S = {1, . . . , n}, /i < st v can be verified in 0(n 2 ) 
time by checking whether [iK < vK holds coordinatewise, where K is the 
n-by-n lower triangular matrix such that K(i,j) = l(i > j) [1]. 

3 Markov processes 
3.1 Markovian couplings 

All state spaces in the following shall be assumed finite or countably infi- 
nite. The keep the presentation short, all results shall be formulated for 
continuous-time Markov processes, which without further mention shall be 
assumed nonexplosive. For a Markov process X with values in S we de- 
note by X(x,t) the value of the process at time t given that was started 
at state x. A Markov process X = {Xi^X^) taking values in S\ x S2 is 
called a Markovian coupling of X\ and X2 if X(x,t) couples X\(x\,t) and 
X2(x2,t) for all t and all x = (x±,X2)- A common approach for showing 
that the time-dependent distributions of two Markov processes X\ and X2 
are stochastically ordered, is to find a Markovian coupling X = (X\, X2) of 
X\ and X2 such that 



almost surely for all t [12]. Observe that if (8) holds, then for all t and all 
increasing functions /, 



0(max(n' 1 , 





(8) 



Ef(X 1 (x u t))=Ef(X 1 (x 1 ,t)) 

<Ef(X 2 (x 2 ,t)) = Ef(X 2 (x 2 ,t)) 
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whenever x\ < x^- Hence Xi(x\,t) < st X2(x2,t). If both X\ and X2 
have unique stationary distributions, it follows by taking limits that the 
stationary distributions are stochastically ordered. 

Observe that (8) is not necessary for the stochastic ordering of the sta- 
tionary distributions of X\ and X2. To formulate a less stringent sufficient 
condition, we shall use the following definitions. Let R be a relation between 
Si and S*2. A pair of Markov processes X\ in S± and X2 in S2 is said to 
stochastically preserve the relation R, if 

xi ~ x 2 => X±(xi,t) ~ st X 2 (x2,t) for all i > 0. 

Moreover, a set B is called invariant for a Markov process X if x G S 
implies G -B for all t almost surely. The following results was proved 

in [9], which we rephrase here for convenience. 

Theorem 3. [9] The following are equivalent: 

(i) X\ and X2 stochastically preserve the relation R. 

(ii) There exists a Markovian coupling of X\ and X2 for which R is in- 
variant. 

Assume now that X\ and X2 are Markov processes with values in an 
ordered space S, and assume ii is a subrelation of the order that is stochas- 
tically preserved by X\ and X2. Then X\{x\, t) < s t X2{x2, t) for all x± ~ X2- 
Especially, if X\ and X2 have unique stationary distributions, then a suffi- 
cient condition for the stochastic ordering of the stationary distributions is 
that X\ and X2 stochastically preserve some nontrivial subrelation of the 
order relation. 

3.2 Subrelation algorithm 

Recall that a matrix Q with entries Q(x, y), x, y G S is called a rate matrix, 
if Q(x, y) > for all x 7^ y and Q(x, x) = — J2 y =/=x Q( x i u) f° r au x - If X is a 
Markov process with rate matrix Q, then Q(x, y) is the transition rate of X 
from state x into y, and we denote by q{x) = —Q(x,x) the total transition 
rate of X out of state x. 

Given Markov processes X\ and X2 with rate matrices Q\ and Q2, define 
the relation-to-relation mapping Mq 1) q 2 by 

M QuQ2 (R) = {(x,y) G R : (v x ,y,v x , y ) G R st }, (9) 
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where the probability measures fj, XtV and v XzV are defined by 

(Jv,y(u) = q~*Qi(x,u) +S x (u), (10) 
v x, v (v) = q-JyQ 2 (y,v) +5 y (v), (11) 

and where q XiV = 1 + qi(x) + 92(2/)- When there is no risk of confusion, we 
denote Mq 1} q 2 = M. Moreover, define recursively the sequence M k (R) by 
setting M°(R) = R, M k (R) = M(M k ~ l (R)) for k > 1, and denote the limit 
of the sequence by 

M*{R) = r\%L M k (R). 

Theorem 4. [9] The relation M*{R) is the maximal subrelation of R that 
is stochastically preserved by X\ and X 2 . Especially: 

(i) Xi and X 2 stochastically preserve R if and only if M(R) = R. 

(ii) X\ and X 2 stochastically preserve a nontrivial subrelation of R if and 
only ifM*(R) / 0. 



Algorithm 2 Computation of R' = M(R). 
R' <— n\-hy-n 2 zero matrix 
for (x, y) € R do 

?<- l + qW(x) +q<?\y) 
fi^q- 1 Q^(x,-)+5 x (-) 

^q- 1 Q^(y,-) + 5 y (-) 

Check whether fj, ~ st v (use Algorithm 1) 

if n ~ s t v then 

R'(x,y)^l 
end if 
end for 



When the state spaces Si and S 2 are finite, Algorithm 2 describes how 
to numerically compute M(R), and Algorithm 3 describes the computation 
of M*(R). Observe that for finite state spaces, Algorithm 3 computes the 
apparently infinite intersection n^L M fc (it!) in finite time, because as long 
as M k (R) and M k ~ 1 (R) are not equal, they differ by at least one element, 
and the sequence M k (R) is decreasing. 



8 



Algorithm 3 Computation of R* = n™ =0 M k (R). 
R' <- M(R) (use Algorithm 2) 
while R' ^ R do 

i?' <- M(i?) (use Algorithm 2) 
end while 



4 Truncation approach 

4.1 Truncation of Markov processes 

If Q is a rate matrix of a Markov process on a countably infinite space S, 
and S n is a finite subset of S, we define the truncation of Q into 5 n by 



Qn(x,y) 



Q(x,y), xj^y, x,y e S„, 

Q(x,y), x = y, x G 5 n . ( 12 ) 



We shall later approximate Q by Q n and use the finite subrelation algorithm 
applied to Q n . To understand the approximation error, we need to study 
how the untruncated process may escape the set S n . Given a rate matrix 
Q on a countable space S, we say that an increasing sequence of finite sets 
S n C S is a truncation sequence for Q, if U^qS^ = S, and 

{y:Q(x,y)>0}c5 n+ i (13) 

for all x E S n . 

Example 2. Let Q be the rate matrix of a Markov process on Z + that is 
skip-free to the right, so that Q(i,j) = for all j > i + 1. Then 5 n = 
Z + n [0, n] is truncation sequence for Q. 

The next result shows that truncation sequences can be constructed for 
most Markov processes encountered in applications. We say that a Markov 
process X with rate matrix Q has locally bounded jumps, if the set {y : 
Q(x,y) > 0} is finite for all x. 

Lemma 1. Any rate matrix Q of a Markov process with locally bounded 
jumps possesses a truncating sequence. 
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Proof. Because S is countable, we may choose an increasing sequence of fi- 
nite sets K n such that U^ =0 K n = S. Using this sequence we may recursively 
define the sets S n by setting So = Kq, and 

S n+ i = S n U K n U J{S n ), n > 0, 

where 

J(S n ) = {y : Q(x, y) > for some x 6 S n } 

denotes the set of states that are reachable from S n by one jump. Then 
U^Q^n = S, because K n C S n+ \ for all n > 0. Moreover, J(S n ) C S n+ ± 
implies that (13) holds for all n, and induction shows that the sets S n are 
finite, because Q has locally bounded jumps. □ 

4.2 Truncation of stochastic relations 

Let R be a relation between countably infinite state spaces Si and S2 ■ If S' t 
is a finite subset of Si, i = 1, 2, we define the truncation of R by 

R' = Rn (S[ x S' 2 ) (14) 

The corresponding stochastic relation R' st , a relation between probability 
measures on the finite spaces S[ and S' 2 , is called the truncation of R s t into 
S[ x S'2. 

If /U, is probability measure on Si having all its mass on S^, we may 
regard probability measure p! { on S[ by identifying subsets of S[ as 

subsets of Si. 

Lemma 2. Let [ii he a probability measure on Si such that ^(S 1 -) = 1, 
i = 1, 2. T/ien (jii, ^2) £ ^st on/y i/ (fj,^, fj,' 2 ) £ R' st . 

Proof. The claim follows using [9, Lemma 5.2], after observing that R' is 
the relation induced from R by the pair (c^i,^), where 0j is the natural 
embedding of S^ into Si , and y! i = \Xi o . □ 

4.3 Truncated subrelation algorithm 

Given a pair of Markov processes Xi and X2 with rate matrices Qi and Q2, 
and a relation i? between Si and S'2, Algorithm 3 describes how to recursively 
calculate R* = Mq q 2 {R) as the limit of the sequence R k = MKq 2 (R). If 
the state spaces are infinite, R k cannot be computed using finite time and 
memory. Nevertheless, when the processes have locally bounded jumps, the 
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truncations T]\[(R k ) of R k into suitable truncation sets Si, at x &2,n can be 
computed precisely, as shall be shown next. 

For any relation R between S± in and S^m ^ M n (R) = Mq 1 ni Q 2 n (i?) 
be the relation given by Algorithm 2 applied to R using the truncations 
Qi )T i and Q2,n of Qi and Q2 as defined in (13). Moreover, denote by M m = 
Mq 1j q 2 the corresponding untruncated mapping. 

Lemma 3. For any pair of Markov processes with locally bounded jumps, 
the truncation of M OQ (R) to S\ tTl x S2 n satisfies 

TniM^R)) = T n {M n+l {T n+l {R))) for all n. 

Proof. Observe that F n (M 00 (i?)) equals the set of points (x, y) € T n (R) such 
that the measures fj, x>y and i/ x>y defined in (10) and (11) are stochastically 
related with respect to R. Because <Si jn and S^n are truncating sequences 
for Qi and Q2, we know that fj, XtV and v x>y have supports in Si >n +i and 
£2,71+1 j respectively. Hence by Lemma 2, T n (M 00 (i?)) equals the set of 
points (x,y) £ T n (R) such that /j, Xty and v X;V are stochastically related with 
respect to T n+ i(R). Because for (x,y) S T n (R), the measures pL X}V and v X) y 
remain the same, if we replace Q\ and Q2 by Qi >n +i and Q2,n+i in (10) 
and (11), the claim follows. □ 

Theorem 5. For any pair of Markov processes with locally bounded jumps, 
the truncation of R k = M^{R) to Si jTl x S2, n satisfies 

T n (R k ) = T n (M n+1 T n+1 ) ■ ■ ■ (M n+k T n+k )(R). 

Proof. Apply Lemma 3 and induction. □ 

Algorithm 4 describes how Theorem 5 can be used to compute T n (R k ) 
in finite time and memory. 

The following result is a necessary condition for finding a subrelation of 
R that is stochastically preserved by a pair of Markov processes. 

Theorem 6. Let X\ and X2 be Markov processes with locally bounded 
jumps. IfT n (M k (R)) = for some k, then T n (M*(R)) = 0. Especially, 
if for all n, T n {M k (R)) = for some k, then X\ and X2 do not stochasti- 
cally preserve any nontrivial subrelation of R. 

Proof. Because the sequence M k (R) is decreasing, and because the trunca- 
tion map is monotone with respect to set inclusion, the first claim follows. 
For the second claim, observe that if T n (R*) = for all n, then because 
U„(Si, n X S 2 , n ) = S, R* = 0. □ 
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Algorithm 4 Computation of R' = T N {M K {R)). 
R! <- T N+K {R) 
for k = 1, . . . , K do 

n <- iV + if + 1 - k 

Qi,n <— truncation of Q\ into Si jTl 

Q2,n truncation of Q2 into S2,n 

R''^T n {R') 

R' <— Algorithm 2 applied to (Qx t n->Q%mR') 
end for 



5 Applications 



5.1 Multilayer loss networks 
5.1.1 Overflow routing 

Consider a loss system with K customer classes and two layers of servers, 
where layer 1 contains m^ servers dedicated to class k, and layer 2 con- 
sists of n servers capable of serving all customer classes. Arriving class-A; 
customers are routed to vacant servers in one of the layers, with prefer- 
ence given to layer 1; or rejected otherwise. For analytical tractability, we 
assume that the interarrival times and the service requirements of class-A; 
customers are exponentially distributed with parameters A^ and respec- 
tively, and that all these random variables across all customer classes are 
independent. For brevity, we denote m = (mi, . . . , mx), A = (Ai, . . . , \r), 
and n = (/xi, . . .,hk)- 



n 



mi 



m 2 



Ai 



m 3 



Layer 2 



Layer 1 



A 3 



Figure 1: Two-layer loss network with three customer classes (K = 3). 
Denote by Xi^(t) the number of class-A; customers being served at layer 
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i at time t. The system is described by the Markov process X = (Xi k) 
taking values in 

Si = {x € if x if : Xl , k < m k Vfc, £f =1 ^2,fc < ™} , (15) 
and having the transitions 

' x + e ljk , at rate X k l(x ljk < m k ), 

x + e 2k , at rate Afcl(xi fe = m fc , ^fe-i £2,fc < n), 
x - ei^, at rate fJ, k x± ik , 
£ - e 2 ,fc, at rate fJ. k x 2>k , 

where ej k denotes the unit vector in Zf x Zf corresponding to the coordi- 
nate direction (i,k). 



5.1.2 Maximum packing 

To approximate the original two-layer loss system, we consider a modifica- 
tion of the system, where customers are redirected from layer 2 to layer 1 as 
soon as servers become vacant. This corresponds to the so-called maximum 
packing policy introduced by Everitt and Macfadyen [4]. Denote by Y k (t) 
the total number of customers in the system with maximum packing. Then 
t \— > Y(t) = (Yi(t), . . . ,Yx(t)) is a Markov process (see [5]) with values in 

5*2 = [y G ■ Ylk=iiVk ~ m k)+ < ™} , 
and having the transitions 

(y + e k , at rate X k l(y + e k € S 2 ), 

y l— > s 

[y-e k , at rate fi k y k . 

The structure of the above transition rates implies that the stationary dis- 
tribution 7ry of Y is a product of Poisson distributions truncated to S 2 [7] , 
so that 

, , A (w^) yfc c o 

k=i Vk - 

where the constant c can be solved from ^ y eS 2 7ry (^) = ^" This product 
form structure allows for fast computation of stationary performance char- 
acteristics of the maximum packing system. 
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5.1.3 Stochastic comparison 

Table 1 illustrates the outcomes of the subrelation algorithm (computed 
using [8]) applied to various initial relations, where 

tf sum = {(x, y)eS 1 xS 2 : £ fc (xi, fc + x 2 , k ) < £ fc Vk} , 
RT m = {(x, y)eS 1 xS 2 : £ fe x 1>k < J2 k (y k A m k )} , 
^coord = y )eSixS 2 : x 1)k < y k A m k for all k} . 

Note that i?| um relates the state pairs (x, y) where the total number of 
layer-1 customers corresponding to x is less than that corresponding to y. 
Moreover, i?™° rd relates the state pairs where the populations in layer 1 are 
coordinatewise ordered. The entry "several" in Table 1 refers to running 
Algorithm 2 separately for several pseudorandom parameter combinations, 
and taking the intersection of the produced relations as a final result. 



R 


(Ai,A 2 ) 




M*(R) 


R snm 
R snm 
R sum 


(1,1) 
several 
several 


(1,1) 

(1,1) 
several 


R sum n ^sum R ^rmn 
j^sum p| ^coord 




R snm 
R snm 


(1,1) 
several 
several 


(1,1) 

(1,1) 
several 


^coord 
^coord 


^>coord 


several 


several 


^coord 



Table 1: Outcomes of the subrelation algorithm for a network with mi = 
1, m 2 = 1, n = 2. 

From Table 1, we can make several conclusions on the behavior of the 
subrelation algorithm: 

• The pair (X\,X 2 ) appears to stochastically preserve i?| um for all pa- 
rameter combinations (same parameters in both systems). This fact 
is in fact not hard to verify analytically. 

• When /ii = /i2, the pair (X\, X 2 ) stochastically preserves a nontrivial 
subrelation of R snm . The maximal subrelation of R sura preserved may 
depend on the model parameters. 

• R sum n Rf ord appears to be a subrelation of R sum that is stochastically 
preserved for all parameter choices, as long as fi\ = [i 2 (Figure 3). This 
fact is proved in [5]. 
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• When the system is symmetric (Ai = A2, fJ>i = ^2 and m\ = rn-2), a 
larger relation R snm n i?f m n Rf in D i? sum n i2f ord is stochastically 
preserved (Figure 2). 

• When fi\ 7^/42) -R sum in general does not have a nontrivial subrelation 
preserved by the pair. This fact is also reflected in [5, Example 5.2.1], 
where it was found that the stationary distributions are not in general 
stochastically ordered with respect to R sum . 

As an illustration, the limiting relations R sum n i?f m n Rf in and i? sum n 
^coord ^fji} ec [ c i rc les), together with the initial relation R snm (filled + unfilled 
circles) are plotter in Figures 2 and 3, respectively. 
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Figure 2: Relations R snm and M*(R sum ) for a system with mi = l,rri2 = 
1, n = 2, for A = (1, 1) and /j, = (1, 1). 



5.2 Parallel queues 

Consider a system of two queues in parallel, where customers of queue k 
have arrival rate A^ and service rate Assuming that all interarrival and 
service times are independent and exponential, the queue length process 
X = (X\, X2) is a Markov process in with transitions 

{x + efc, at rate A&, 
x — e/c, at rate /ifcl(xfc > 0). 

We shall also consider a modification of the system, where load is balanced 
by routing incoming traffic to the shortest queue, modeled as a Markov 
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Figure 3: Relations R snm and M*(R sum ) for a system with mi = l,rrt2 = 
1, n = 2, for various A and fixed fi = (1, 1). 



process X LB = (X^ 6 ,^ 5 ) in I? + with transitions 

at rate (Ai + A2)l(^i < X2) + Ail(xi = X2), 
at rate (Ai + X 2 )l(xi > x 2 ) + A 2 l(3;i = x 2 ), 
at rate nil(x\ > 0), 
at rate H2l(x 2 > 0). 

Common sense suggests that load balancing decreases the number of 
customers in the system in some sense. The validity of this comparison 
property can be numerically studied using Algorithm 4. Denote the rate 
matrix of X LB by Q± and the rate matrix of X by Q 2 , and let S n = 1?+ n 
[0,n — l] 2 , n > 1. Then S n is a truncation sequence for both Q\ and Q 2 . 

Figure 4 illustrates five iterations of the subrelation algorithm (computed 
using [8]) truncated to S% applied to the coordinatewise order 

flcoord = | (z>y) e z 2 + xZl: Xl < yi , x 2 < y 2 } , 

with the parameters Ai, X 2 , 1^1,^2 all equal to one. Because T^(R 4 ) = 
T3(i? 5 ) = 0, we conjecture that there exists no nontrivial subrelation of 
i? coord stochastically preserved by (X hB ,X) (see Theorem 6). 
Let us next study another order on Z^_, defined by 

R snm = {(x,y) eZ 2 + xZ 2 + :x 1 +y 1 , x 2 + y 2 } . 

Figure 5 illustrates five iterations of the subrelation algorithm (computed 
using [8]) truncated to S3 applied to R° = R sum with the parameters 



x 



' x + ei, 

X + 32, 

x - ei, 
x - e 2 , 
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Figure 4: Truncated subrelation algorithm applied to R 



coord 



Ai, A2, Hi, fJ>2 an equal to one. The observation that T^(R k ) remains un- 
changed from k = 1 onwards suggests that some nontrivial subrelation of 
R sum might be stochastically preserved by (X LB ,X). Indeed, it has been 
analytically shown [9] that whenever fii = M2, the untruncated subrelation 
algorithm converges to the relation 



R* 



{(x, y) : xi + x 2 < yi + y 2 and x\ V x 2 < yi V y 2 } . 



As a consequence of Theorem 4, the pair (X LB ,X) stochastically preserves 
the relation R*, which may be identified as the weak majorization order on 
1% [10]. Especially, 

X\ B {t) + X^ B (t) < st Xi(t)+X 2 (t), 
xl B (t) V X 2 LB (t) < st Xi{t) V x 2 (t), 

for all t, whenever the initial states X LB (0) and X(0) satisfy the same in- 
equalities. 
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Figure 5: Truncated subrelation algorithm applied to R s 



6 Conclusions 



This paper presented computational methods for verifying stochastic re- 
lations and finding relation-invariant couplings of continuous-time Markov 
processes on finite and countably infinite state spaces. A key point of the 
paper is that the stochastic relationship between two probability measures 
can be quickly numerically checked, if one of the measures has small support 
(Theorem 2). This result allows the development of a truncation approach 
for finding relations stochastically preserved by pairs of Markov processes 
with locally bounded jumps. The truncated subrelation algorithm (Algo- 
rithm 4) allows to numerically find candidates for a subrelation of a given 
relation that is stochastically preserved by a pair of Markov processes. It 
remains an interesting open problem for future research to study how the 
truncated subrelation algorithm behaves for structured Markov processes 
with for example shift-invariant transition rate matrices. 
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